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INTRODUCTION 



The latest tally of multiple plane t candidate systems dis- 
covered by the Kepler mission (iBatalha et alj bOld) m 



ABSTRACT 

A fraction of multiple planet candidate systems discovered from transits by the Kepler 
mission contain pairs of planet candidates that are in orbital resonance or are spaced 
slightly too far apart to be in resonance. We focus here on the four planet system, KOI 
730, that has planet periods satisfying the ratios 8:6:4:3. By numerically integrating 
four planets initially in this resonant configuration in proximity to an initially exterior 
cold planetesimal disk, we find that of the order of a Mars mass of planet-orbit-crossing 
planetesimals is sufficient to pull this system out of resonance. Approximately one 
Earth mass of planet-orbit-crossing planetesimals increases the interplanetary spacings 
sufficiently to resemble the multiple planet candidate Kepler systems that lie just 
outside of resonance. This suggests that the closely spaced multiple planet Kepler 
systems, host only low mass debris disks or their debris disks have been extremely 
stable. We find that the planetary inclinations increase as a function of the mass in 
planetesimals that have crossed the orbits of the planets. If systems are left at zero 
inclination and in resonant chains after depletion of the gas disk then we would expect 
a correlation between distance to resonance and mutual planetary inclinations. This 
may make it possible to differentiate between dynamical mechanisms that account for 
the fraction of multiple planet systems just outside of resonance. 

iLibert fc Tsiganid I2OIII : iRein et al.l I2OI2I ). Gravitational 
interactions between planets leading to planet-planet 
scattering events (e.g., ^Gozdzic wski fc Mi gaszcwski' 200£; 



eludes 361 multiple-planet svstems /Pab rvckv et al.f2012h . 
A statistical analysis, focused on the probability that bi- 
nary stars are the most likely contaminant, finds that most 
of the multiple plane t candidates are real planetary systems 

■ (|Lissauer et al.ll2012l '). The large number of recently discov- 
. 5^ ] ered multiple planet systems represents a significant (by an 

■ order of magnitude) increase in the number of known multi- 
5^ ■ pie planet syste ms compared to thos e discovered from radial 
Ci , velocity surveys IWright et al.ll201ll ). 

Both transit and radial velocity discovered multiple 
planet systems contain p airs of planets that are in first order 
mean motion resonance I Wright et al.ll2009l : II issauer et al.l 
I2OIII : IWright et al] l201ll 'l. In the Kepler transit systems, 
there are statistically significant excesses of candidate planet 
pairs both in resonance and space d slightly too far ap art to 
be in resonance (as delineated by IVeras fc FordI 20121'). par- 



Fabrvckv fc Murrav-Clavi .2010. ; iMoore fc Quillen ,20121: 



Rein et al.l |2012| ). turbulence in the disk (jPierens et al 



2011 
2011 



and tidal interactions between planets (jPapaloizou 



ticularly near the 2:1 mean motion resonance ( Lissauer et al.l 
l201ll : lFabrvckv et al.|[20ll ). 

Planet migration due to tidal interaction with a gas 
disk is a possible mechanism through which convergent 
migration ind uces resonance capture, leaving planets in 
resonance ( e.g. Lee fc Peald 2002 ^: Fcrraz-McUo ct al. '2003'; 
Kiev et all |2004|: iLed 12004 iPapal oizou fc Szuszkicwicz 



20051 : iThommes et al.l I2OO8I : fMorbidelli et" al] feoOTl : 



Lithwick fc Wull2012l ) have been proposed as mecha- 
nisms for pulling initially resonant planetary systems out 
of resonance. Interactions with planetesimals, for example 
as part of the 'Nice' model for the early solar system 
evolution, can also cause plan ets to diverge awa y from 
resonance (|Tsiganis et al.l bOOSl ) (also see iThommes et al.l 
l2008h . 

Simulations of planets and planetary embryos em- 
bedded in a g as disk allow pla n ets to become t r apped 
in resonance (iKlev et al.l |2004| : iMorbidelli et al.l l2007l : 
I Rein et "all |2012| V After the gas disk dissipates, newly 
formed pla nets may be left in a chain of mean-motion 
resonances (IMorbidelli et al.l l2007l : iMatsumura et al.l I2OI0I : 
iMoeckel fc Armitagd I2OI2I '). By a chain of mean motion 
resonances, we mean that each consecutive pair of plan- 
ets is in a J -I- 1 : j, first order, mean motion resonance 
(though the integer j can differ for each pair). Resonant 
chains have been chosen as initial conditions for studies 
of planetary system evolutio n afte r the depletion of the 
gas disk (e.g., Tsiganis e t al] l2005l : IThommes et all l2008l : 
iBatvgin fc Brown! 1201^ When pairs of planets are in or 
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near mean motion resonances, there may be a librating or 
nearly fixed Laplace angle involving three or more bodies. 
Three-bod y resonances ( e.g., those comprised of zero-th or- 
der terms; iQuillenllioTll ) may be present even when pairs 
of planets are not in first order mean motion resonances. A 
system in a resonant chain of first order resonances does not 
necessarily exhibit a librating Laplace angle involving three 
or more bodies. 

The four planet candidate system, KOI-730, contains 
planets with periods that satisfy the ratios 8:6:4:3 to ap - 
proximately 1 part in 1000 or better (|Lissauer et al.|[201ll ). 
These ratios imply that the outer two and inner two planet 
pairs are in (or near) a 4:3 mean motion resonance and that 
the middle pair of planets is in (or near) a 3:2 resonance. In 
this study we explore the evolution of a similar model four 
planet system in proximity to a planetesimal disk. We ask: 
how much mass in orbit- crossing planetesimals is sufficient 
to pull this system out of resonance? 

For a pair of p lanets and a first order mean motion res- 
onance, [Pabrxcky^et al. (2012) define a parameter, ("i.i, that 
measures proximity to a first order mean motion resonance. 



Table 1. Masses and Periods for the planet candidates 
in the KOI 730 system and initial orbital elements 



Round 



(1) 



where V = Pi/Pj (greater than 1) is the observed period 
ratio of planets ^, ,?' (and as defined in the appendix by 
iFabrvckv et al.ll2012l 'l. The 'Round' function rounds to the 
nearest integer. The parameter C,i.i = at a j -|- 1 : j first 
order resonance and the function goes from -1 to 1 at sec- 
ond order resonance s {j -\- 2 : j). The parameter used by 
iLissauer et al.l ()201ll ). ("i = l-5Ci,i and varies from -1.5 to 
1.5. 

The probability density distribution of ("i values gen- 
erated from all pairs of Kepler planet transit candidates, 
for pairs residing in a single system, exhibits a peak at 
about Ci f« —0. 2, (see Figure 11 by iLissauer et all I2OIII 
and Figure 5 bv lFabrvckv et al.ll2012l ). For KOI-730, Ci = 
-0.0123, -0.0186, -0.0063 for the inner p air, middle pair 
and outer pair of planets, respectively (|Fabrvckv et al.l 
|2012| ) , consequently the system is likely to be in or very near 
a resonant chain. By integrating a system modeled after the 
KOI-730 system that is initially in a chain of resonances, we 
ask: how much mass in orbit-crossing planetesimals would 
increase the interplanetary spacing so that i^i ~ —0.2? We 
also note that due relatively low mass of the planets and 
their proximity to the star, their resonant widths are likely 
smaller than ("1 = 0.1, as first order mean motion resonance 
width sc a les ap proximately with mass to the 2/3 power (e.g. 
IWisdomI (|l980l ')'). 

We first describe how we find resonant chain configura- 
tions for the KOI-730 system. We use these configurations 
to construct initial conditions for N-body integrations that 
include a planetesimal disk. A summary and discussion fol- 
lows. 



2 SETTING UP RESONANT CHAIN 
CONFIGURATIONS 

We first describe how we find orbital elements consistent 
with a 8:6:4:3 ratio resonant chain for a four planet system 
similar to the KOI-730 multiple planet system. We find res- 



Mass 


Period 


a 


e 




M 


(A%) 


(days) 










2.5 


7.3840 


1.0 


0.055589 


-2.808830 


1.216544 


3.7 


9.8487 


1.211221 


0.071128 


-0.739956 


2.841707 


8.6 


14.7884 


1.586171 


0.050110 


1.453201 


-0.924252 


6.2 


19.7213 


1.920630 


0.043428 


-1.990222 


-0.172501 



Planet masses for the KOI 730 system are computed 
from radii based on transit durations taken from 
http: / / archive. stsci.edu /kepler /planet_candidat es.html| 
l|Batalha et al.l l20lj ~ using equation |3] The observed peri- 
ods are given in days and are taken from the same website. 
The rightmost four columns give the orbital elements we used 
as initial conditions (see section 2). Here ui is the argument of 
pericenter and M the mean anomali. Initial inclinations and the 
longitudes of the ascending node were set to zero. These orbital 
elements were taken from the integration shown in Figure [T] 



Table 2. N-body Simulations 



Simulation 

Z 

M 

E3 

E 

N5 

N 



Planetesimal Mass 

(Me) 


10-* 
3.3 X 10-4 

10"^ 
3.3 X 10-3 
1.67 X 10-2 



Orbit-crossing mass 


0.04 
0.12 
0.46 
1.7 
16.6 



Each simulated disk contains 1024 equal mass planetesimals. The 
planetesimal masses in units of Earth mass are listed in the second 
column. The third column shows the total mass in planetesimals 
(in Earth masses) that crossed the orbits of any of the planets 
at the end of the simulation. The names of the simulations are 
related to the masses of the planetesimal disks. The Z simulation 
has a zero mass disk. The M, E and N simulations have disks 
with Mars, Earth and Neptune masses, respectively. The E3 and 
N5 simulations have disks with a third Earth and a fifth Neptune 
mass, respectively. 

onant chain configurations by integrating a 5 body system 
under the influence of gravity (four planets and the central 
star) and including a Stokes drag-like form for dissipation 
that induces both migration and eccentricity damping (as 
previously done by 'Beauge et al.l I2OO6I : iBatygin fc Brownl 
[2OIO; Libe rt fc Tsig anis 2011]). T he drag gives a forc e per 
unit mass in the form adopted bv lBeauge et al.l (|2006[ ) 

F - - ^-^'^ (')\ 

r drag — 9 

where v is the planet velocity and Vc is the velocity of a 
planet in a circular orbit at the current radius (from the star) 
of the planet. We use a 4th o rder adaptive step-size He rmite 
integrator (that described bv lMakino fc Aarseth|[l993 ) with 
the addition of the above drag force. We work in units of 
the innermost planet's initial orbital period or 7.3840 days. 

The drag force induces radial migration on a timescale 
I ~ Ta, and eccentricity damping on a timescale | ~ 
where a,e are the planet's semi-major axis and eccentric- 
ity, respectively. Resonant capture can only occur when the 
drift rate, a, is sufficiently slow such that the square of li- 
bration frequency in resonance exceeds the drift rate (this 
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defines the adiabatic limit; iQuillenlliooi ). FoUowine reso- 
nant capture, eccentricities of pla nets can increase as they 
drift inwards. jLee fc Pealell2002l). causing instability (e.g., 
iKlev et all [20041 : iLibert fc Tsiganisll2"oill : lRein et"alll2012h . 
For two planet systems, an equilibriu m state may be rea ched 
that depend s on the ra t io K = ^ (|Lee fc Pealdl2002l V As 
discussed bv lRein et al.l ()2012l ). if two bodies lie initially out- 
side the 2:1 or 3:2 resonance, then they are more likely to 
capture in one of those than in the 4:3 resonance. Here we 
require that the outer and inner pair of planets are cap- 
tured into the 4:3 resonance, consequently we began the 
integration with planets spacings just outside the desired 
resonances. The migration rates were adjusted so that the 
migration is sufficiently fast that capture into weaker second 
order resonances such as the 5:3 or the 7:5 is unlikely. 

The masses of the four planets were set from the 
radii measured from the transit durations and reported by 
iBatalha et ai] (|2012l ). We adopt the power-law relationship 
for planetary mass, Mp, as a function o f radi us, Rp, used by 
iLissauer et al] (l201ll ): iFabrvckv et all l|2012D . 



(3) 



where , Rg) are the mass and radius of the Earth, and 
the exponent a — 2.06 for Rp > R^. This choice is moti- 
vated by Solar System planets as it is a good fit to Earth, 
Uranus, Neptune, and Saturn. KOI-730's surf ace gravity and 
effect ive temperature (reported in Table 9 bv lBatalha et al] 
|2012| . with surface gravity logiO g(in cgs) = 4.39 and effec- 
tive temperature Tsff = 5599 K) are similar to that of the 
Sun (with logio=4.43 and Te// = 5780K) so we computed 
planet to stellar mass ratios from the estimated planet radii, 
equation [3] and using a Solar mass for the host star. Esti- 
mated planet masses and their periods are summarized in 
Table [T] The estimated planet masses for the KOI-730 sys- 
tem are much lower than the approxim ately Jupiter mass 
planets considered bv lRein et al] (|2012l l for the HD2006964 
system. The 4:3 resonance may be more stable in lower mass 
systems. 

An integration is shown in Figure[l]with a) showing the 
four semi-major axes as a function of time and b) showing 
the period ratios of consecutive pairs of planets. Semi- major 
axes are shown in units of the initial innermost planet's semi- 
major axis. The timescales Ta and can be chosen sepa- 
rately for each planet. We set Ta for the innermost planet 
to be extremely long (10^ orbits), and that for the outer 3 
planets to be progressively shorter ranging from 10^ to 10''' 
orbits. Note that 10^ orbits of the innermost planet is only 
approximately 2000 years. We arranged the drift rates so 
that the outer planet migrates more quickly than the inner 
ones. After the outer planet captures the third planet, the 
two together migrate more slowly than the outer planet. To 
maintain a constant drift rate for the pair, the third planet 
was set with a longer value for Ta, and similarly, Ta for the 
second planet was chosen to be larger than that for the third 
planet. 

Once in resonance, the eccentricities of the planets in- 
crease as the planets drift inwards. A steady state can be 
reached with eccentricity values that depe nd on the size of 
the eccentricity damping, set here with Te (|Leell200^ ). High 
values of Te (corresponding to low levels of damping) are 
associated with high planet eccentricities, whereas low val- 
ues of Te reduce the eccentricities of the planets. We found 




40000 



80000 120000 
time (orbits) 



160000 



Figure 1. a)Semi-major axes of the four planets under the in- 
fluence of a forced dissipative process as a function of time, b) 
Period ratios for the same integration. This integration was used 
to generate initial orbital elements for the four planets in the sim- 
ulations including planetesimals. We took orbital elements from 
the time at 70, 000 orbits. 



that this system became unstable without significant eccen- 
tricity damping. Consequently we set the level of eccentric- 
ity damping sufficiently high, Te ^ 10'', so that the system 
remained in resonance after all planets were captured into 
resonance. This value is high compared to that predicted for 
tidal interactions between disk and planet. Large K values 
(up t o 100) have been used previously (Bat y gin fc B r ownl 
|2010D to generate stable initial resonant conditions for sub- 
sequent integration, as we do here. As can be seen from 
Figure [T] the outer pair of planets first captures into reso- 
nance, then the middle pa ir and lastly the innerm ost pair 
is captured into resonance. iMorbidelU eral] (|2007| ) stressed 
that the order of captures can affect the resonant angles. 
The 8:6:4:3 ratios imply that the second and fourth planets 
are near or in a 2:1 mean motion resonance and the first and 
third planets are also near such a resonance. We find that 
fine-tuning in initial planet semi-major axes, forced drift and 
eccentricity damping rates are required to put the system in 
the 8:6:4:3 chain of resonances. 

We are not investigating formation mechanisms (see 
I Rein et al ]|2012l for a first investigation into this tricky prob- 
lem). 



3 N-BODY SIMULATIONS WITH A 
PLANETESIMAL DISK 

Using orbital elements from the converging integration dis- 
cussed above, we run N-body simulations of the four planets, 
initially in resonance, in the vicinity of a planetesimal disk 
and about a central star. These simulations were run with 
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the software QYMSYM (jMoore fc Quillenll201ll 'l. QYMSYM 
is a GPU-accelerated hybrid second order symplectic N- 
body integrator which permits close encounters similar to 
the M ERCURY software package developed bv IChamberj 
l|l999t ). 

We began each simulation with four planets in a reso- 
nant chain and an external planetesimal disk comprised of 
1024 equal mass planetesimal particles. Orbital elements for 
the four planets in the integration shown in Figure [1] at time 
t = 70, 000 orbits, and listed in Table [U were used as initial 
conditions for these integrations. 

The planetesimal disk particles ranged in semi-major 
axis from amin ~ 1.95, just outside the outermost planet, to 
dmax = 2.95. The distribution of planetesimal semi-major 
axes is flat with probability independent of semi-major axis 
within Umin and Umax- The initial eccentricity and inclina- 
tion distributions were chosen using Rayleigh distributions 
with the mean eccentricity e equivalent to twice the mean 
value of the inclination i and i = 0.01. The initial orbital 
angles (mean anomalies, longitudes of pericenter and lon- 
gitudes of the ascending node) for the planetesimals were 
randomly chosen. We work in units of the initial orbital pe- 
riod of the innermost planet. 

Six simulations were run, each with different planetesi- 
mal mass. We included a test case with massless disk parti- 
cles to check the stability of the integration lacking planetes- 
imals. We labeled the simulations by the mass of the plan- 
etesimal disk, with simulation Z corresponding to a massless 
planetesimal disk while simulations M, E, and N correspond 
to Mars, Earth and Neptune mass disks. Simulation E3 and 
N5 refer to simulations with a third of an Earth mass and 
a fifth of a Neptune mass disk, respectively. Each simula- 
tion was run for 500,000 orbital periods (or 10^ years as the 
innermost planet's orbital period is only 7.4 days). 

The timestep was chosen to be 0.08607 out of a possible 
2n orbit. Given this step size, energy conservation (measured 
by dE/E) was 10~^ or better for all simulations except dur- 
ing the simulation labeled N (its energy conservation was 
8 X 10"^). 

Due to the proximity of the planets of KOI- 730 to their 
star, the planets practically fill their Hill sphere. The large 
fraction of the Hill sphere filled limits the escape velocities 
of close encounters. Our integrator does not check for col- 
lisions though it does integrate close encounters. To take 
into account the large fraction of the Hill sphere filled, we 
adjusted the smoothing lengths during close encounters be- 
tween planets and planetesimals. 

In these simulations we used a smoothing length of s = 
1 X 10""^, which corresponds to a distance that is slightly 
larger than the radii of the innermost planet (about 1.1 x). 
The Hill radius of the innermost planet, which would be the 
smallest of the four due to it both being the closest to the 
central star as well as the least massive, is a little more than 
ten times larger (about 13.5 x) than this smoothing length. 

We have checked that the inclination distributions and 
extent of migration are not strongly dependent on the as- 
sumed smoothing length. This was done by comparing the 
results of our simulation described above that have planet 
radii sized smoothing lengths to a separate set of identical 
simulations only with a smoothing length 100 times smaller. 

During each simulation we computed the mass in plan- 
etesimals that crossed the orbits of the planets. We count 



a planetesimal as orbit-crossing if its pericenter is less than 
the apocenter of the outermost planet. The total planetesi- 
mal mass that crossed the orbits of the planets at the end of 
the simulations are also listed in Table [21 As we will discuss 
below, these masses can be used to place limits on the total 
quantity of planetesimals that may have crossed the orbits 
of planets in the KOI- 730 system. 

3.1 Planetary migration in the N-body 
simulations 

Figure [2^) shows period ratios for each consecutive pair of 
planets as a function of time for all 6 N-body simulations. 
For the inner and outer planet pairs (initially in 4:3 res- 
onance) the period ratio is shown subtracted by 4/3. The 
middle planet pair (initially in 3:2 resonance) is plotted sub- 
tracted by 3/2. Our test case Z simulation (bottom panel in 
Figure [2^) remains stable throughout the integration, con- 
sequently we are confident that the orbital elements chosen 
from our capture integration are stable. It is possible that 
this system becomes unstable on a timescale longer than 
500,000 orbits. 

The M, E3, E, N5, and N simulations contain disks with 
increasing mass. Because of the proximity of the outermost 
planet to the inner edge of the planetesimal disk, the outer 
planet migrates outwards as planetesimals cross the orbits 
of the planets and exchange angular momentum with them. 
The outer planet migrates furthest in the simulations with 
highest disk mass. For the most part, period ratios increase, 
though as the planets separate, mean motion resonances, as 
they are crossed, can cause jumps in both eccentricity and 
semi-major axis. At some times we see a signature of three- 
body resonances, when the period ratio for one consecutive 
pair (of three planets) increases as the period ratio of the 
other consecutive pair decreases l|Quillenll201ll ). 

The Z simulation remains locked in resonance through- 
out the integration. However, we see that the M and E3 
simulations no longer maintain their chains of MMR's at 
250k and 175k orbits respectively. The E simulation is re- 
moved from resonance at approximately 50k orbits. This 
allows us to estimate the mass in planet-crossing planetesi- 
mals required to pull this system out of resonance for each 
simulation, finding on average rric ~ MMara- The N and N5 
simulations have their planets out of resonance very quickly 
due to the larger amount of orbit crossing mass early in the 
simulation. In both of those simulations the planets are out 
of resonance too quickly for us to accurately determine how 
much mass was orbit crossing. Simulations E, E3, and M all 
have total crossing masses within an order of magnitude of 
each other. 

For each integration we compute C,\^\ (equation [l]) as a 
function of time for each consecutive pair of planets. Figure 
shows the function computed from all consecutive 
planet pairs for the same simulations. Variations in the 
parameter are largest for the most massive disk. We see from 
Figure [2l3 that C,i.\ approaches -0.2 at a time of about 50k 
orbits for planets 3 and 4. This allows us to estimate the 
amount of planetesimal disk mass that would move a sys- 
tem sufficiently out of resonance to contribute to the peak 
seen in the C,i distribution of the Kepler multiple planet can- 
didates. Based on the mass in orbit- crossing planetesimals 
in simulation N5 we estimate rric « is required to pull 
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a system initially in a resonant chain sufficiently far apart to 
give a (i equivale nt to the position of th e peak seen in the 
(^1 distribution by iFabrvckv et al.l (|2012h (see their Figure 
5). 

We note in Figure [^Jj that we find that the (i^i function 
does not always decrease for each planet pair. Furthermore, 
during migration we see a range of ("1,1 values. To produce 
a peak in the distribution of ("1 we would require a very 
specific or fine tuned value for the mass in orbit- crossing 
planetesimals for multiple planet systems that are originally 
in resonant chains. 

Figure [2j; shows the inclinations of all planets as a func- 
tion of time. We see that the planet inclinations increase to a 
couple of degrees in the simulations with the most massive 
disks. We have checked that simulated planetesimals have 
not been ejected at high velocity. The masses of our simu- 
lated planetesimals is quite low (see Table [T]) and the incli- 
nations slowly increase. The increases in planet inclination 
are unlikely to be numerically generated during encounters 
and due to extreme scattering events. The smooth increase 
in inclinations can be attributed to a combination of gravi- 
tational heating caused by scattering of planetesimals and to 
crossing of vertical resona nces resulting from the migration 
of the planets (as seen bv lLibert fc TsiganidlioTll '). 

Figure [2^ implies that planet inclinations in closely 
spaced multiple planet systems depend on the total plan- 
etes imal mass that has cro ssed the orbits of the planets. 

iFabrvckv etlH (|2012h find that mutual inclinations for 
the multiple planet systems lie in the range 1.0 — 2.3° and 
that their distribution is well modeled by a Rayleigh dis- 
tribution with standard deviation ^ {P) = 1.8°. We see in 
Figure[2j; that the inclinations rise above this value when the 
total mass in orbit crossing planetesimals exceeds 17M® , or 
about one Neptune mass. Therefore, closely packed systems 
similar to the KOl-730 system likely have not experienced 
an era similar to the Late-Heavy Bombardment in our Solar 
system. From this comparison we tentatively place a limit 
on the total mass in orbit crossing planetesimals during the 
lifetime of the KOI-730 system of rric < M^. 



4 DISCUSSION AND CONCLUSION 

Using an N-body integrator with the addition of drag caus- 
ing convergent migration and eccentricity damping, we con- 
structed a resonant chain for four planets with period ratios 
8:6:4:3, and used planet masses estimated for the KOI-730 
multiple planet system. Orbital elements from this integra- 
tion were then used as initial conditions for N-body sim- 
ulations with the four planets which included an external 
planetesimal disk. Interactions with the planetesimal disk 
allowed the planets to migrate, primarily diverging rather 
than converging, as expected. 

We find that one Earth mass of orbit-crossing planetes- 
imals is sufficient to pull a system similar to the KOI-730 
system out of its chain of mean motion resonances. As planet 
radii estimated from transit data depend on stellar radii, it 
is possible that we have over or underestimated the masses 
of the planets in the KOI-730 system. The distance migrated 
by a planet should scale with the mass in planetesimals that 
it interacts with. However, mean motion resonant widths are 
larger when planet masses are larger. We find that it is more 




SOOOO UXKXK) 150000 200000 250000 300000 350000 400000 450000 500000 
orbital periods 




50000 100000 150000 200000 250000 300000 350000 400000 450000 500000 
orbital periods 



50000 100000 150000 200000 250000 300000 350000 400000 450000 500000 
orbital periods 

Figure 2. a) Period ratios for consecutive pairs of planets sub- 
tracted by their initial value and as a function of time. Each panel 
shows period ratios for a different simulation. From top to bottom, 
the planetesimal disk mass decreases. Red, green and blue lines 
refer to the inner, middle and outer consecutive planet pairs, re- 
spectively, b) (^11 parameters computed from consecutive planet 
period ratios and for the same simulations. The colors are the 
same pairs as in a) . c) Planet inclinations in degrees as a function 
of time for the same simulations. The red, green, blue and pink 
lines refer to the first (inner), second, third and fourth (outer) 
planets, respectively. 



© 0000 RAS, MNRAS 000, 000-000 



6 



difficult to form a resonant chain via resonant capture and 
sucli chains are less stable when the planets are more mas- 
sive. Consequently, the amount of material required to pull 
a system out of resonance may not scale linearly with planet 
meiss. Nevertheless, we expect that if the true planet masses 
are larger than adopted here, a larger mass in orbit- crossing 
planetesimals would be required to pull this system out of 
resonance and vice- versa if the planets are less massive. 

After the system is out of resonance, we find 
that it remains stable. Our simulations do not exhibit 
planet/planet orbit crossing events. The KOI- 730 system, 
comprised of sub-Neptunian mass planets, can be com- 
pared to the HR8799 system, comprised of hyper- Jovian 
mass planets in resonance. When the HR8799 system is 
'ulled out of resonance, the system is extremely unstabl e 
Gozdziewski fc Migasze^skH l2009l : iMoore fc Quillenll2012l ) . 

The KOI-730 system is likely in (or very close to) a 
resonant chain. Because interactions with planetesimals and 
tidal interactions between planets primarily cause planetary 
orbit s to slowly diverge (.Papaloizou 2011; Lithwick & Wu| 
|2012| ) and so pull systems out of resonance, we can be fairly 
confident that processes following depletion of a gas disk 
did not put this system in resonance. Because the system is 
currently near or in resonance, only a small mass in plan- 
etesimals could have ever crossed the orbits of these planets. 
We infer that this system either lacks a debris disk or con- 
tains one that is so diffuse or stable that less than an Earth 
mass of debris has ever crossed the orbits of these planets. 

We also find that an Earth mass of orbit-crossing plan- 
etesimals can cause the planets to migrate far enough that 
the system lies sufficiently outside of resonance to resem- 
ble the Kepler system s with Ci —0.2 where the peak 
of th e distribution lies (|Fabrvckv et al.ll2012l : iLissauer et al.l 
I2OIII ). However we expect that different planetary systems 
would have different quantities of orbit-crossing planetesi- 
mals. Consequently we would not expect that a narrow peak 
in the C, distribution would arise in a distribution of plane- 
tary systems. Fine tuning in the quantity of planetesimals 
may be required to account for the peak seen in the ( distri- 
bution. The tidal dampin g scenario for pulling pairs of plan- 
ets a way from resonance (jLithwick fc "Wull2012l : [Papai oizoul 
I2011I ) may more naturally account for a peak in ^. 

We find that the more an initially flat planetary system 
interacts with planetesimals, the higher the mean planet 
inclinations. If somewhere between a few earth masses to 
about a Neptune mass of planetesimals cross the orbits of 
the planets, the planet inclinations can increase to a few 
degrees. This is sufficiently high that they would not be all 
simultaneously be detected in transit. The tidal circulariza- 
tion scenario would not give a relation between migration 
distance and inclination. However over short distances, as 
planets migrate and they cross vertical resonances, we would 
expect a correlation between migration distance (and so dis- 
tance out of resonance) and inclination. Study of the relation 
between the inclination and period distributions of the Ke- 
pler systems may differentiate between roles of tidal forces 
and planetesimals. 

It is possible that the transiting multiple planet sys- 
tems discovered by the Kepler mission are more compact or 
lower mass than radial velocity discovered planetary sys- 
tems. Can we differentiate between the tidal interaction 
mechanism for pulling systems out of resonance and that 



caused by interactions with planetesimals? If planetesimals 
interact with planets, then it is likely that planet inclina- 
tions can increase as planets cross vertical resonances. It 
may be possible to differentiate between these two mecha- 
nisms based on inclination distributions as tidal interactions 
likely do not increase inclinations but planet /planetesimal 
interactions can. Note that Libert & Tsiganis (2011,) have 
shown that resonant capture for higher planet mass systems 
can also induce inclination variations. However, additional 
mechanisms, such as turbulence associated with a gas disk or 
secular perturbations from distant planets could also affect 
the inclination distributions. Future studies can probe the 
role of planetesimals, and migration associated with scatter- 
ing them, in accounting for the inclination distributions of 
the Kepler planetary systems. 

We have focused here on interactions with a low ec- 
centricity planetesimal disk. When it encounters a planet, 
a high eccentricity planetesimal is less strongly gravitation- 
ally focused than a low eccentricity one. Consequently, high 
eccentricity objects are less effective at scattering a planet 
or inducing migration. Our limit on the total mass in or- 
bital crossing planetesimals can be considered a lower limit 
as we began with a low eccentricity disk. Future studies can 
explore the possibility that compact Kepler systems har- 
bor massive outer planetary systems and high eccentricity 
cometary populations. 

In summary, we believe that closely spaced, low incli- 
nation multiple planet Kepler systems likely have either low 
mass or extremely stable debris disks. There appears to be 
a relation between the inclination and amount of migration 
for a planet. The inclination distributions may make it pos- 
sible to differentiate between dynamical scenarios for pulling 
planets out of resonance. Due to the improbability of a Late 
Heavy Bombardment like scenario for KOI-730, we believe 
these inclinations are probably caused by crossing vertical 
resonances. 
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